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We present a detailed analysis of binary black hole evolutions in the last orbit, and demonstrate 
consistent and convergent results for the trajectories of the individual bodies. The gauge choice can 
significantly affect the overall accuracy of the evolution. It is possible to reconcile certain gauge 
dependent discrepancies by examining the convergence limit. Wc illustrate these results using an 
initial data set recently evolved by Briigmann et al. ( Phys. Rev. Lett. 92, 211101). For our highest 
resolution and most accurate gauge, we estimate the duration of this data set’s last orbit to be 
approximately 59A/adm- 
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Introduction. Over the course of the next decade, 
instruments capable of detecting gravitational radiation 
(such as LIGO, VIRGO, TAMA, GEO600) are expected 
to open a new observational window on the universe. The 
collision of binary compact objects such as black holes 
(BHs) is one of the most promising sources for first gen¬ 
eration gravitational wave observatories. The theoretical 
framework for modelling binary BH systems is the com¬ 
plete set of nonlinear Einstein equations. Intensive efforts 
to develop numerical codes able to solve these equations 
using supercomputers, have shown that it is now possible 
to evolve BHs for periods of an orbit Ull If these 
simulations are to produce waveforms useful for detector 
searches, high demands are placed on their accuracy [j. 

The near-zone dynamics of binary BH systems are no¬ 
toriously difficult to simulate, and to analyse. Using the 
BSSN formulation and a particular set of gauges, a series 
of BBH configurations, corresponding to initial data in 
quasi-circular orbit at successively larger separations @, 
were all found to coalesce in slightly more than a half 
orbit 1- A similar BSSN evolution carried out using 
somewhat different gauges and numerical methods for 
another data set, slightly further out along the orbital 
sequence, was found to evolve for much more than the 
estimated orbital timescale 114Madm without finding a 
common apparent horizon (AH) |l|. In fact, no common 
horizon was found long after the BHs would reasonably 
be expected to have merged. 

In this paper, we carry out an evolution of the same 
data set and show that it does indeed carry out a com¬ 
plete orbit before a common AH forms. As the BH sepa¬ 
ration decreases, a local measure of the angular velocity Ul 
increases, so that the duration of the final orbit is approx¬ 
imately 59M. The trajectories are convergent for a range 
of resolutions and within a class of gauge conditions. 

However, we do find that very high resolutions are re¬ 
quired in order to obtain evolutions close to the con¬ 
tinuum limit. The resolutions we have applied here are 
significantly higher than those used in analogous BH evo¬ 


lutions to date, except for Q where similar resolutions 
were used. With insufficient resolution, we show that 
it is possible to substantially under- or over-predict the 
orbital period. We also find that apparently small devia¬ 
tions in the chosen coordinate conditions (gauge), based 
on choice of parameters within a particular family, can 
have a strong influence on the discretisation error within 
the subsequent evolution. 

Methods. Initial data for the evolutions discussed 
in this paper correspond to Brandt-Briigmann “punc¬ 
tures” [6]. The particular orbital parameters are chosen 
to be identical to those first evolved in Q, namely initial 
separation L/M = 9.32 M, bare masses m = 0.47656M 
of each BH, and equal and opposite linear momenta, 
p = ±0.13808. These parameters are chosen to approx¬ 
imate a BH pair in quasi-circular orbit. The angular 
velocity H 0 = 0.054988 suggests an orbital timescale of 
114M for a perfectly circular orbit. The initial geome¬ 
try is determined by numerically solving the constraint 
equations using the solver of Ansorg et al. 

The binary BH evolution is carried out using the 
“BSSN” formulation of the Einstein equations @, B, Eoj, 
with an implementation described explicitly in llj. We 
use a version of the “l±log” lapse and T-driver shifts, 
given by Eqs. (33) (with f(a) = 2A>^ L /a) and (46) (with 
F(cx) = 3/4a p /'I>B L ) of fn). We generalise the parame¬ 
ter rj to a lapse-dependent coefficient a q r), where typically 
q € [0,4]. This is helpful for long-term evolution of BHs 
from large initial separations through merger, and will 
be discussed in detail elsewhere. 

We also dynamically adapt our gauges to the hori¬ 
zon location so as to approximate co-motion. Individ¬ 
ual horizons are located often during the evolution using 
the finder described in 0- From the motion of the hori¬ 
zon centroid, the angular and radial positions and veloci¬ 
ties are determined. Using a damped harmonic oscillator 
equation, the acceleration needed to bring the horizon 
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centroid back to the initial position is determined by 

A = -[2TQA+(A-A°)]/T 2 . (1) 

Here A stands for either the current azimuthal angle <fi 
or radius r = x 2 + y 2 in the orbital plane, and A 0 for 
its initial value. T and Q are constants determining the 
timescale and damping factor of the harmonic oscillator. 
The correction is added to the time derivative of the shift 
vector as 

A/3* = {-y, x, Or 6 - {x, y, 0) 1 r/r°. (2) 

The punctures are excised using an extension of the 
“simple excision” techniques which have proven success¬ 
ful in evolving single BH spacetimes |u| [l4| . In partic¬ 
ular, we apply the boundary condition to an embedded 
boundary whose shape is determined by the horizon lo¬ 
cation, as described in mu 

Spatial differentiation is performed via straightfor¬ 
ward finite-differencing, incorporating nested mesh- 
refined grids with the highest resolution concentrated in 
the neighbourhood of the individual horizons. The mesh 
refinement is implemented via the Carpet driver Jld for 
Cactus. The evolutions carried out in this paper made 
use of 8 levels of fixed 2:1 refinement. We fix the re¬ 
gions of increased resolution around the initial BH lo¬ 
cations. Because we are making use of a BH-adapted 
gauge, the individual horizons remain on the fine grids 
throughout the evolution without requiring moving grids 
or excised regions. We have used finest grid resolutions of 
ft = 0.025 M, 0.02M, 0.018M, 0.015M and 0.0125M, with 
an outer boundary at 96M in all cases. Our overall finite 
differencing accuracy is second order in space and time. 

As a diagnostic of the dynamics of the individual BHs, 
we measure the proper distance within a slice between the 
pair of apparent horizons (AHs). This is calculated by 
shooting spacelike geodesics from the origin (taking ad¬ 
vantage of the spacetime symmetry) to one of the horizon 
surfaces 0 As this measure of distance is within each 
slice, it does depend on the particular lapse condition (as 
will be seen below), but still provides useful information 
about the near-zone BH motion in the given slicing. 

Results. Evolutions were carried out for a number of 
resolutions and gauge parameters. A first observation, 
important in interpreting the results of Hi , is that 
the infall coordinate trajectory at a given resolution de¬ 
pends crucially on the gauge choice. In Fig. the re¬ 
sults of evolutions of three choices of gauge parameters 
(introduced above) are displayed. The first, which we la¬ 
bel “GC1”, sets (m,n,p,q,r],T,Q) = (0, 2,1,1,4, 5,1). 
The second gauge choice, “GC2”, sets parameters to 
(4, 2,1,4, 2, 5,1). The third gauge choice, “GC3”, is de¬ 
scribed below. Plotted are the proper separations of the 
AHs within the slice versus the coordinate time. 

We note that for the GG1 evolutions at a grid spacing 
of ft. = 0.025M, the BHs fall together rather rapidly. By 
coordinate time t = 75 M, the proper distance between 
the AHs is down to L/M = 4.5. Increasing the resolution 


to ft = 0.020 M, ft = 0.018M, and ft = 0.015M there is a 
trend towards longer evolution times before reaching the 
same separation, though the convergence is slow. For the 
highest resolution which we have evolved, ft = 0.0125, the 
infall to the same separation is delayed to only t = 88 M. 

The evolutions resulting from gauge choice GC 2 have 
an entirely different character. At the lowest resolution, 
with ft = 0.025M, we find that not only does the system 
evolve well beyond an orbital time period before a com¬ 
mon AH appears, the orbit appears to be elliptical, with 
proper separation first increasing with time, and then 
falling back, reaching a separation L/M = 4.5 only af¬ 
ter t = 119 M. This behaviour is qualitatively similar to 
what was reported for the same initial data set in jl] . In¬ 
creasing the resolution as above, however, we find a faster 
coalescence as resolution is increased, and the apparent 
elliptical nature of the orbit disappears. 

Importantly, the two families of evolutions correspond¬ 
ing to GG1 and GC 2 do show a tendency towards 
each other as resolution is increased. In fact, a 3-term 
Richardson extrapolation (eliminating 2nd and 3rd or¬ 
der error terms) of the computed trajectories show that 
the two families converge to a common result in the con¬ 
tinuum limit. Fig. |T] also shows the results of the ex¬ 
trapolation for both sets of evolutions at the 3 highest 
resolutions where we had sufficient data available. (Us¬ 
ing other sets of 3 resolutions gives very similar results.) 

The differences between GG1 and GC 2 results may at 
first seem to be the expected result of the variation in 
slicing condition between the two classes of evolutions. 
However, the fact that the results converge to the same 
coordinate trajectory in the continuum limit suggests a 
somewhat different explanation, namely that the finite 
difference error inherent in the evolution (particularly in 
the lapse) can be strongly influenced by the gauge choice. 

The Richardson extrapolation provides an indication 
of the accuracy of the evolutions at a given resolution, or 
conversely, the resolution required to obtain a result of a 
given accuracy. The Richardson-extrapolated coordinate 
trajectory seems robust, but we note that the highest 
resolution (h=0.0125M) GC 2 simulations which we have 
carried out still represent a 22% error in separation at t = 
100M compared to the extrapolated result. The sensitive 
dependence of accuracy on gauge suggests that “ideal” 
gauge conditions may be able to improve the accuracy 
greatly at a given resolution. 

Based on results of the previous experiments, a third 
gauge choice parameters was found to demonstrate this 
point. Under the label GG3, evolutions of the same data 
were carried out using gauge parameters (4,2,1,4,2,5,2.6), 
i.e. the same as for GC 2 except for the drift-correction 
damping parameter. The GC 3 results are also plotted in 
Fig. 1. Once again, we find the proper separation mea¬ 
sure to be convergent with increasing resolution, and fur¬ 
ther that the Richardson extrapolation lies on top of the 
curves predicted by the previous evolutions. However, in 
this case the highest resolution attempted, ft. = 0.015M, 
lies much closer to the limiting trajectory. 
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FIG. 1: The minimal proper distance between individual 
AHs as a function of time. Lines show two representative 
resolutions, h = 0.025 M and h = 0.015M for each gauge 
choice GG1, GC 2, and GC 3 discussed in the text, and an 
h = 0.0125 evolution for GCT. Points show Richardson ex¬ 
trapolations (RE) for resolutions h G {0.018, 0.015, 0.0125}M 
for GG1, h G {0.020,0.018,0.015}M for GG2, and 
h G {0.025, 0.020, 0.015}M for GG3. 
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FIG. 2: Angular motion of individual AHs as a function of 
time (lower scale) and angle (upper scale, in radians measured 
backwards along the orbit from the first appearance of a com¬ 
mon AH). The curve shows the angular velocity 12 of the AH 
centroid for GG3, as measured by integrating the corrections 
applied to the gauge via 0. The horizontal line corresponds 
to the fi 0 = 0.055 estimated for the initial data. 


Fig. El plots the evolution of the individual BHs’ an¬ 
gular velocity, as measured by the AH centroid’s drift, 
for GC 3 at our highest evolved resolution. From our ex¬ 
perience with more closely separated binaries, the initial 
shift profile corresponded to a quasi-rigid rotation with 
H = 0.3, approximately 60% of the nominal Ho = 0.055 of 
the initial data 0 . The co-rotation shift adjustment 0 
quickly raises the effective angular velocity of the co¬ 
ordinate system up to a value close to the initial data 
prediction Ho- As expected, this angular velocity then 
gradually increases as the BHs spiral together. 

Fig .□plots the overall motions of the AHs for the same 
GC3 case as fig. 0 The individual horizon shapes are 
shown at intervals of 5M, transformed according to the 
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FIG. 3: Schematic showing the motion of one of the BHs 
with time, for the GG3 gauge choice at the highest resolution 
h = 0.015M. At intervals of t = 5 M, the AH cross-sections 
in the xy plane are plotted by transforming the co-rotating 
coordinate system by the specified angle and distance. The 
apparent growth of the AHs with time is a non-physical co¬ 
ordinate effect. The first appearance of a common AH at 
t = 124M, and corresponding final single AH, are shown su¬ 
perposed on the figure as dotted lines. 


radial and angular motion determined above. The initial 
rapid dip in separation from 9.32 M to 8.3 M in Fig.0can 
now be seen to be largely a result of the initial coordinate 
expansion of the individual horizons, due to our choice 
of zero radial shift in the initial data. 

A common AH is first detected at t = 124AL, by which 
time an angular displacement of 10.2 radians, or 1.6 rev¬ 
olutions, has taken place. The common horizon is found 
via the method of pre-tracking, in which a family of sur¬ 
faces with smallest possible generalised expansion is fol¬ 
lowed to provide an estimate which converges on the first 
common AH 0. The listed time is expected to corre¬ 
spond to the first genuine appearance of a common AH, 
independent of the search algorithm. Counting back¬ 
wards from the appearance of the common horizon, the 
duration of the last orbit is approximately 59M. It is 
interesting to note that for the cases where a common 
horizon is found, it occurs when the proper separation is 
approximately 1.8M. Our Richardson-extrapolated GC1 
and GC 2 evolutions reach this same separation less than 
5M earlier than this GC 3 evolution. 

Given the demonstrated resolution dependence of these 
results, it is likely that these results will be subject to 
some modification as more accurate evolutions become 
available. It should also be recalled that a merged event 
horizon within the slice will have formed earlier than the 
common AH. Experience with closely separated binaries 
suggests this is typically not more than about 5 M before 
the appearance of the first common AH, but studies of 
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the event horizon evolution for these spacetimes will need 
to be the subject of future studies. We also note that 
the notion of an “orbit” for closely separated BHs is an 
intrinsically gauge-dependent quantity - for example, a 
sufficiently small lapse in the region of the horizons could 
be used to indefinitely delay merger. The slicings used 
here are quite similar in profile to maximal slices, and 
thus not atypical for numerical relativity simulations. 

As a final point, we note that a number of measures 
of accuracy have been monitored during the course of 
these evolutions. In particular, the constraint violation 
is found to remain below a value of 0.05 at all points 
outside the BH horizons and away from the boundaries 
for the duration of the runs. The AH masses of the in¬ 
dividual BHs were measured and found to be essentially 
constant, convergent and accurate even for the low reso¬ 
lution runs. We have also compared binaries with closer 
separation, such as the “QC-0” data evolved by a num¬ 
ber of groups and find excellent agreement 

with published results. 

Conclusions. We have carried out evolutions of BBH 
configurations, from initial data in quasi-circular orbit, 
through plunge, to the formation of a common horizon. 
The dynamics exhibit a number of interesting properties. 

Using a measure of angular velocity based on the AH 
motion, we have found that the black holes remain sep¬ 
arated for more than 1.5 revolutions. This measure ini¬ 
tially yields an angular velocity which is consistent with 
the initial data model. The most accurate estimate of 
the duration of the last complete orbit before formation 
of the first common AH was 59 M. As might be expected 
for bodies falling towards each other, this orbital period 
is considerably less than is predicted by the initial an¬ 
gular velocity, 114M. The results have been carried out 


at resolutions much higher than similar studies to date, 
and exhibit good numerical consistency under a range of 
resolutions and for a variety of gauge parameter choices. 

By Richardson-extrapolating the black hole trajecto¬ 
ries we derived what appears to be a robust extimate 
of their continuum limit, which proved to be almost in¬ 
dependent of gauge choice, at least within the family of 
gauges considered here. 

These evolutions required extremely high resolution to 
attain good accuracy. Insufficient resolution can result in 
very different predictions for the orbital trajectories and 
the period of the final orbit. This period directly influ¬ 
ences the phase of asymptotically observed wave forms 
at their strongest point, and thus is crucial to reproduce 
accurately. 

We demonstrated that the gauge choice can feed di¬ 
rectly into the numerical accuracy of the solution, due 
to unwanted dynamics in the evolution variables caused 
by gauge effects. A preferred gauge was found to sig¬ 
nificantly improve the accuracy at a given resolution. A 
more detailed exposition of the gauge conditions used 
here, applied to a sequence of initial data configurations, 
will be the subject of upcoming studies. 
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